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Abstract 

For compressible fluids under shock wave reaction, we have proposed two Multiple-Relaxation- 
Time (MRT) Lattice Boltzmann (LB) models [F. Chen, et al, EPL 90 (2010) 54003; Phys. Lett. A 
375 (2011) 2129.]. In this paper, we construct a new MRT Lattice Boltzmann model which is not 
only for the shocked compressible fluids, but also for the unshocked compressible fluids. To make 
the model work for unshocked compressible fluids, a key step is to modify the collision operators of 
energy flux so that the viscous coefficient in momentum equation is consistent with that in energy 
equation even in the unshocked system. The unnecessity of the modification for systems under 
strong shock is analyzed. The model is validated by some well-known benchmark tests, including 
(i) thermal Couette flow, (ii) Riemann problem, (iii) Richtmyer-Meshkov instability. The first 
system is unshocked and the latter two are shocked. In all the three systems, the Prandtl numbers 
effects are checked. Satisfying agreements are obtained between new model results and analytical 
ones or other numerical results. 
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I. INTRODUCTION 



In recent years, the Lattice Boltzmann (LB) method has attracted much attention as a 
powerful tool in direct numerical simulation of fluid flows l|-|3| . Unlike traditional computa- 
tional fluid dynamics methods which solve macroscopic governing equations, the LB method 
employs the discrete Boltzmann equation which describes the fluid on the mesoscale level. 
This kinetic nature provides the LB method with essential physics. 

However, there are also some limitations that restrict the applications of traditional LB 
method, such as the numerical stability problem, the fixed Prandtl number, and so on. The 



stability problem has beenpartly addressed 



method 



r pro l 



flux-limiterp] and dissipation 



jy a number of techniques, such as the entropic 
% |8[ techniques. Besides these techniques, an 



effective method is the Multiple Relaxation Time (MRT) LB method which employs 

multiple relaxation parameters in the collision step, instead of the commonly used Single 
Relaxation Time (SRT) collision. The flexibility gained from the MRT collision can be used 
to improve the stability property and overcome the fixed Prandtl number problem. 



Tot 

system 



re authors' knowledge, most of the existing MRT LB models work only for isothermal 
15|, to cite but a few. To simulate system with temperature field, Luo, et al-flf]] 



12 



suggested a hybrid thermal MRT LB model, in which the mass and momentum equations 
are solved by the MRT model, whereas the diffusion-advection equation for the temperature 
is solved by Finite Difference (FD) technique or other means. Guo, et al.[l^ proposed a 
coupling MRT LB model for thermal flows with viscous heat dissipation and compression 



work. Mezrhab, et al.[18[ proposed a double MRT LB method, where MRT-D2Q9 model and 
the MRT-D2Q5 model are used to solve the flow and the temperature fields, respectively. 
Besides the models mentioned above, we have proposed two MRT finite difference Lattice 



Boltzmann models for compressible fluids under shock in previous work[19|, |20(. Numerical 
experiments showed that compressible flows with strong shocks can be well simulated by 
these models. In this paper, we further propose a new MRT Lattice Boltzmann model, 
which is not only for the shocked compressible fluids, but also for the unshocked compressible 
fluids. The rest of the paper is organized as follows: In Sec. II, we present the MRT LB 
model. The von Neumann stability analysis is given in Section III. Simulation results are 
presented and analyzed in Section IV. Section V makes the conclusion. 
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FIG. 1: Schematics of Vj for the discrete velocity model. 



II. DESCRIPTION OF THE MRT LB MODEL 

In the MRT LB method, the evolution of the distribution function /j is governed by the 
following equation 

^+^ = -M^,(A-/f), (1) 

where Vi a is the discrete particle velocity, i — 1,. . . ,N, N is the number of discrete velocities, 
the subscript a indicates x or y. The variable t is time, x a is the spatial coordinate. The 
matrix S = MSM -1 = diag(si, s%, •• • , s/v) is the diagonal relaxation matrix, /j and /j 
are the particle distribution function in the velocity space and the kinetic moment space 
respectively, f\ = mijfj, rriij is an element of the transformation matrix M. Obviously, the 
mapping between moment space and velocity space is defined by the linear transformation 
M, i.e., f = Mf, f = M _1 f, where the bold-face symbols denote N-dimensional column 
vectors, e.g., f = (/i,/ 2 ,--- Jn) T , f = (A, A,-" Jn) T , M = (m u m 2 ,--- ,m N ) T , = 
(run, m>i2, • • • , mix). f^ q is the equilibrium value of the moment /j. 

We construct a two-dimensional MRT LB model based on a 16-discrete- velocity model 
(see Fig. 1): 

'eye: (±1,0), forl<i<4, 
(±1,±1), for5<i<8, 
eye : (±2, 0) , for 9 < % < 12, 
(±2, ±2) , for 13 < i < 16, 
where eye indicates the cyclic permutation. 



{vii,v a ) 
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The transformation matrix M is constructed according to the irreducible representation 
bases of SO (2) group, and it can be expressed as follows: 



where 



M = (mi,m 2 , 



m x = (1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1), 
m 2 = (1, 0, -1, 0, 1, -1, -1, 1, 2, 0, -2, 0, 2, -2, -2, 2), 
m 3 = (0, 1, 0, -1, 1, 1, -1, -1, 0, 2, 0, -2, 2, 2, -2, -2), 
"4 = (^,^,1,1,1,1,2,2,2,2,4,4,4,4), 
m 5 = (1, -1, 1, -1, 0, 0, 0, 0, 4, -4, 4, -4, 0, 0, 0, 0), 
m 6 = (0, 0, 0, 0, 1, -1, 1, -1, 0, 0, 0, 0, 4, -4, 4, -4), 
m 7 = {\, 0, -i 0, 1, -1, -1, 1, 4, 0, -4, - 

m 8 = (0,i 0,-^,1,1, -1,-1,0, 4, 0,-4 



, ^, vj, "i u/| 



m 9 



1, 0, -1, 0, -2, 2, 2, -2, 8, 0, -8, 0, -16, 16, 16, -16), 



m w = (0, -1, 0, 1, 2, 2, -2, -2, 0, -8, 0, 8, 16, 16, -16, -16), 

mi1 = 4' i' i' I' 4 ' 4 ' 4 ' 4 ' 16 ' 16 ' 16 ' 16 )' 

m 12 = (1, 1, 1, 1, -4, -4, -4, -4, 16, 16, 16, 16, -64, -64, -64, -64), 

m 13 = (1, -1, 1, -1, 0, 0, 0, 0, 16, -16, 16, -16, 0, 0, 0, 0), 

Tom = (0, 0, 0, 0, 2, -2, 2, -2, 0, 0, 0, 0, 32, -32, 32, -32), 

m 15 = (1, 0, -1, 0, -4, 4, 4, -4, 32, 0, -32, 0, -128, 128, 128, -128), 

m 16 = (0, -1, 0, 1, 4, 4, -4, -4, 0, -32, 0, 32, 128, 128, -128, -128). 

For two-dimensional compressible models, we have four conserved moments, density p, 
momentums j x , j y , and energy e. They are denoted by fx, f 2 , f$ and / 4 , respectively. 



Specifically 



expansion 
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h = jx, h = jy, U = e = P(T + u 2 /2). Using the Chapman-Enskog 
2l| on the two sides of LB equation, the Navier-Stokes (NS) equations 



for compressible fluids can be derived. The equilibria of the nonconserved moments can be 
chosen as 

./r U; --.!;):,>■ (2a) 
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^ = J X Jy/ P , (2b) 

f? = (e + pRT) Jx /p, (2c) 

/f = (e + P RT)j y /p, (2d) 

if = VI - (2e) 

f$ = (V 2 x -jl)j v /f?, (2f) 

f? 1 =2f/p-Vl+jlW, (2g) 

Al = (6pe - 2j* - 2#(£ - j, 2 )/ P 3 , (2h) 

f3 = (6pe-2g-2fi)j x j y /ft>. (2i) 

The recovered NS equations are as follows: 

^ + ^£ + % = o, (3a) 

f - 1 K/o) - 1 (Wrt ~£ + 1 w£ - If )1 + 1 W£ + > Pb) 

f + 1 («./,) h- I = -f + + £>] - £w£ - |f)i. (30 



l + l [(e+p)wrf+ » [(e+iu/rf 



where /x s = pRT/s 5 , pi v = pRT/s 6 , X 1 = 2pRT/s 7 , A 2 = 2pRT/s 8 . 
When p s = p v = /i, Ai = A2 = A, the above NS equations reduce to 



Mr - ^ 

dj a | dQajp/p) OP | g <% Q | 
<9t eta;/? <9o; Q dx@ dxp dx a dx x a 
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It should be pointed out that, the viscous coefficient in the energy equation ( Hcj) is not 
consistent with that in the momentum equation fflbl) . Motivated by the idea of Guo et al. 

n 

[171 ]. the collision operators of the moments related to the energy flux are modified: 
S 77 (/ 7 - fT) =► MA - f?) + (s 7 /s 5 - l) P Tu x (^ - ^) + (s 7 /s 6 - l)pTu y (^ + 

S 88 (/ 8 - fT) => S 88 (/ 8 - f?) + (s 8 /s 6 - l)pTu x (^ + ^) + (V* - l )P Tu v^ - Tif)- 
With this modification, we are able to get the following t her mohydro dynamic equations: 



(5a) 

Ul ux a 

9ja . 9 Unj.i/p) __ ■ | ,-■'..„ 

~E + dx t —I" 1 — — ''-''I- 

| + J_ |(e + P)jM = _y Afl _^_ + + g _ ^ (5c) 

The 



dt dx a 


= o, 














dP d 

dx a dxp 


OX 13 


dx a 


du x 
dx x 
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This modification method is also suitable for our previous MRT models 
definitions of /if, /if, /if have no effect on macroscopic equations, so the choices of the 
three moments are flexible. Now we give three different and typical formations: /^| = /^f = 
fH = (version 1); = 0, Q = pu x (-4 + l0T + 5u 2 x -5u 2 y ), f% = pu y (4-10T+5u 2 x -5u 2 y ) 
(version 2); /ff = M 12i fr\ ft! = M 15i f^, f% = M 16i fr x , fT* = p/^T) exp(-(^ Q - 
u a ) 2 /(2T)) (version 3). In the second version, the MRT model reduces to the usual lattice 
BGK model in ref. [22] which uses a higher-order velocity expansion for Maxwellian-type 
equilibrium distribution, if all the relaxation parameters are set to be a single relaxation 
frequency s, namely S = si. 



III. STABILITY ANALYSIS 



In this section, the von Neumann stability analysis [20j on the new MRT LB model is 
performed. In the stability analysis, we write the solution of FD LB equation in the form of 
Fourier series. If all the eigenvalues of the coefficient matrix are less than 1, the algorithm 
is stable. Coefficient matrix Gij of the unmodified model can be expressed as follows, 

n — A __ Via ^ (Ak a Ax a _ -ik a Ax a \c , 1 , VioA& s2, ik a Ax a _ 9 

+ e- a ^)8 ij - AtM^S ifc (|| - (6) 
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where 



fk — M fc p/p, 



df? _ dfk Q dp + djfdT + df? du a 



(7) 



dfj dp dfj dT dfj du a dfj ' 
Because of the modification of the collision operators, the coefficient matrix Gij correspond- 
ing to energy flux should be replaced by 



Via At 

2Ax n [ 



2 V Ax f 



,df 7 df 



cq 







du 



77K df 3 

dUr, 



d 



dfj ' dt 



[(s 7 /s 5 - l)pTu x 



, du r du 



dx dy 



V -)] 



dy 



(8) 



and 



G a So 



V in At 



2Ax a 



)dij + 2 [ AxJ [ 



+ e 



ik Q ^A cc qi 



)^-AtMr 8 i{S 88 (-^ 



<■<] 







du,, du 1 



dfj df, ' df 



— [(s 8 /s 6 - l)pT^(-^ + -^)] 



<9x dy 



d On 



du,, 



)]}■ 



(9) 



dx dy 

We conduct a quantitative analysis using the software, Mathematica. In Fig. 2 we show 
some stability comparisons for the new MRT model, its SRT counterpart and our previous 



model in refs. 
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201 ] . The abscissa is for kdx, and the vertical axis is for |co>| ma :z which is the 



largest eigenvalue of coefficient matrix G^. Grid sizes are dx = dy = 10 , and time step 
is dt = 10 -5 , the relaxation frequency in SRT is s = 10 5 . The other parameters in stability 
analysis are chosen as follows: (a), (p, Ui, u 2 , T) = (2.0, 2.0, 0.0, 2.0), the collision parameters 
in MRT are s { = 10 5 ,i = 1, • • • , 16, the Mach number is 1 (Ma = u/VZT = 2/2); (b), 
(p, «i, u 2 , T) = (2.0, 6.0, 0.0, 2.0), the collision parameters in the three versions are s 9 = 10 3 , 



those in model 



191 ] are s 10 = 5 x 10 4 , s n = 2 x 10 4 , Si 3 = 1.5 x 10 4 , and those in |20j] are 



sg = 8 x 10 3 , S13 = 7 x 10 4 , S14 = 5 x 10 4 , the others are 10 5 , and the Mach number is 
3.0; (c), (p, Ui,u 2 , T) = (2.0, 10.0, 0.0, 2.0), the collision parameters in the three versions are 
s 9 = 1.2 x 10 4 , Si 3 = 10 2 , S14 = 5 x 10 4 , si 5 = 1.5 x 10 3 , those in modelfl9| are s 13 = 1.5 x 10 4 , 
and those in model|20] are sg = 2 x 10 3 , Sis = 6.1 x 10 4 , su = si$ = 3 x 10 4 , the others 
are 10 5 , and the Mach number is 5; (d), (p,ui,u 2 ,T) = (2.0,12.0,0.0,2.0), the collision 
parameters in the three versions are s 9 = 10 4 , s ri = 10 2 , s i4 = 6 x 10 4 , s i5 = 1.5 x 10 3 , 




FIG. 2: Stability comparison for the new MRT model and its SRT counterpart. 



in Q 



and those in model 19] are su = 2 x 10 4 , srs = 1.5 x 10 4 , and those in |2(| are s$ = 10 3 , 
S13 = 5 x 10 3 , Si4 = S15 = 3 x 10 4 , the others are 10 5 , and the Mach number is 6. In case 
(a), the MRT and SRT have the same stability; with the increase of Mach number (case (b) 
and (c)), the MRT models are stable, while the SRT version is not; if further increases the 
Mach number, MRT models also encounter instability problem (case (d)). It is clear that, 
by choosing appropriate collision parameters, the stability of MRT can be much better than 
the SRT. Three versions of the new MRT model do not show large differences in numerical 
stability. 

IV. NUMERICAL SIMULATIONS 



In this section, we study the following problems using the modified MRT LB model: 
Couette flow, One-dimensional Riemann problem, and Richtmyer-Meshkov instability. We 
work in a frame where the constant R — 1. 



S 




FIG. 3: Slip velocity and temperature jump simulated with the three versions and the model 
proposed by Kataoka, et al. 



A. Unshocked compressible fluids 



Here we conduct a series of numerical simulations of Couette flow. The aims of simulation 
of Couette flow are twofold. At first, we prove the Maxwellian property of the discrete 
equilibrium functions. Consider a viscous fluid flow between two parallel flat plates, moving 
in the opposite directions, U wr = —U w \ = 0.2, where subscripts wr and wl indicate the walls 
in the right and left sides. The initial state of the fluid is p — 1, T — 1, U — 0. The 
temperatures of walls are T wr = T w i = 1. Near the walls, we adopt the diffuse reflection 
boundary conditions proposed by Sofonea, et al[23[. In the other two boundaries the periodic 
boundary condition is adopted. In the diffuse reflection boundary, the particles leaving the 
wall are assumed to follow the Maxwellian distribution. Following the discretization of the 
velocity space, in the FD LB model the Maxwellian distribution function is replaced by the 
equilibrium distribution function. 

Fig. 3 shows the velocity and temperature profiles simulated with the three versions of this 



proposed model and the model proposed by Kataoka, et al[24j. The abscissa ix is the index 
of lattice node in the x- directions, and the vertical axes are velocity u and temperature 
T, respectively. The parameters are dx = dy = 0.01, dt = 10~ 4 , NX x NY = 100 x 5. 
The diffuse reflection boundary conditions work well with our model, the slip velocity and 
temperature jump near the walls are clearly seen, and increase with Knudsen number. While 
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FIG. 4: Effects of heat conductivity on temperature profiles of Couette flow, (a) corresponds to 
the unmodified model(version 1), and (b) corresponds to the modified model. Pr = 0.01, Pr = 1 
and Pr = 2 correspond to S7 = s$ = 10, sj = s$ = 10 3 , and S7 = ss = 2 x 10 3 , respectively (other 
collision parameters are 10 3 ). 

it fails to work for the model by Kataoka et al, because the temperature near the wall is 
lower than the wall temperature, which is contrary to physical idea. In Couette flow the fluid 
at the walls should have a higher temperature than the walls themselves, because of the heat 
generated by the viscous flow. We think this contradiction is caused from the equilibrium 
distribution function in their model which is not a Taylor expansion of the Maxwellian. 
So it departs from the basic assumption of diffuse reflection boundary. None of the three 
versions violates the basic assumption and destroys the Maxwellian property of the discrete 
equilibrium functions. 

Secondly, we will compare the ability of the unmodified model and the modified model for 
the unshocked compressible fluids. In the simulation, the left wall is fixed and the right wall 
moves at speed U = 0.1. The simulation results are compared with the analytical solution: 

T = T 1 + (T 2 -T 1 )| + ^ 2 |(1-|), 

where 7\ and T 2 are the left and right wall's temperatures (T\ = 1, T 2 = 1.005), H is 
the width of the channel. Other parameters remain unchanged. Periodic boundary con- 
ditions are applied to the bottom and top boundaries, and the left and right walls adopt 
the nonequilibrium extrapolation method. Fig. 4 and Fig. 5 show the temperature profiles 
of Couette flow simulated with the unmodified model (version 1) and its modified version. 
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FIG. 5: Effects of viscosity on temperature profiles of Couette flow, (a) corresponds to the unmod- 
ified model(version 1), and (b) corresponds to the modified model. Pr = 10, Pr = 5 and Pr = 1 
correspond to S5 = sq = 10 2 , S5 = sq = 2 x 10 2 , and S5 = sq = 10 3 , respectively (other collision 
parameters are 10 3 ). 

In Fig. 4, we fix viscosity coefficient S5 = = 10 3 , and change the thermal conductivity 
s 7 — s 8 from 10 to 2 x 10 3 . On the contrary, we fix thermal conductivity Sj = s§ = 10 3 , and 
change the viscosity S5 = Sq from 10 2 to 10 3 in Fig. 5. (a) corresponds to the unmodified 
model (version 1), and (b) corresponds to the modified model. It is clearly shown that the 
simulation results of modified model are in agreement with the analytical solutions, and the 
Prandtl number effects on unshocked compressible fluids are successfully captured by the 
modified model, but not by the unmodified model. 



B. Shocked compressible fluids 

(a) Riemann problem 

Here we construct a high Mach number shock tube problem with the initial condition, 

(p, Ul ,u 2 ,T)\ L = (5.0,45.0,0.0,10.0), x < 0. 
(p,u 1 ,u 2 ,T)\ R = (6.0,-20.0,0.0,5.0), x > 0. 

The Mach number of the left side is 10.1 (Ma = u/V2T = 45/v^0), and the right is 6.3 
(Ma = u/\/2T = 20/vTO). Figure 6 shows the comparison of LB results and exact solutions 
at t — 0.018, where the parameters are dx = dy = 0.003, dt = 10~ 5 , s 5 = s 6 = 1.5 x 10 4 , 
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FIG. 6: LB results and exact solutions for shock tube problem at time t = 0.018. p: density, P: 
pressure, U : the x— component of velocity, T: temperature. 

other values of s are 10 5 . Squares correspond to simulation results with the unmodified 
model (version 1), the circle symbols correspond to the modified MRT simulation results, 
and solid lines represent the exact solutions, respectively. It can be seen that the simulations 
of the two MRT models do not show large differences. For shocked compressible flows, there 
exist a fast procedure and a slow one. The shock dynamic procedure is fast, while that of 
heat conduction is slow. In such a case, from the view of macroscopic description, the terms 
related to viscosity and heat conductivity may be neglected. So, terms related to viscosity 
and heat conductivity in Eqs.(j3]) and (jSJ) are all small terms and make negligible effects. 
That is the reason why the unmodified model works also well in such cases, 
(b) Richtmyer-Meshkov instability 



The Richtmyer-Meshkov (RM) instability 25|, |26] is a fundamental fluid instability that 
develops when an incident shock wave collides with an interface between two fluids with 
different densities. This instability is involved in numerous physical processes, such as in- 
ertial confined fusion, supersonic and hypersonic combustion, supernova explosion, and so 
on. RM instability has attracted considerable attention for several decades because of its 
important theoretical and practical significance. To the best of our knowledge, the research 
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FIG. 7: SnapslfOts of shock wave reaction orf^hrgle bubble. The left column^(a) corresponds to 
s$ = sq = 10 4 , st = ss = 10 5 , the middle column (b) corresponds to s§ = sq = 10 4 , sj = ss = 10 , 
and the right column (c) corresponds to S5 = sq = 10 4 ,S7 = s$ = 10 3 . From black to white the 
grey level corresponds to the increase of density. 

of RM instability by LB method is still very limited. In this paper, we study the thermal 
conductivity and viscosity effects on RM instability with the MRT LB method. 

The investigation of the interaction of a planar shock with an isolated gas bubble is of 
special significance in the study of RM instability, because the interface of gas bubble has 
typical three dimensional characteristic and large initial distortion. It helps to understand 
the mechanism of RM instability process. The problem we simulated is as follows: A planar 
shock wave with the Mach number 1.22 (D = 1.725), traveling from the right side, impinges 
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FIG. 8: Snap^ts of shock wave reaction ngle bubble. The left column c ^a) corresponds to 
s$ = sq = 10 3 , sj = ss = 10 5 , the middle column (b) corresponds to s§ = sq = 10 3 , sj = ss = 10 , 
and the right column (c) corresponds to S5 = sq = 10 3 ,S7 = s$ = 10 3 . From black to white the 
grey level corresponds to the increase of density. 

on a cylindrical bubble. The initial macroscopic quantities are as follows: 



(p,u 1 ,u 2 ,p) \ x ,y,o-- 



(1,0,0,1), pre — shock, 

;i.28, -0.3774, 0,1.6512), post - shock, (11) 
(0.1358,0,0,1), bubble, 

The domain of computation is a rectangle Nx x Ny = 600 x 100, Nx and Ny are the 
numbers of lattice node in the x- and y- directions. Initially, the bubble is at the position 
(450,50), the post-shock domain is [501,600] x [0,100]. In the simulations, the right side 
adopts the initial values of post-shock flow, the extrapolation technique is applied at the left 
boundary, and reflection conditions are imposed on the other two surfaces. 
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In Figures 7 and 8, we show some simulation results with different configurations. The 
abscissa is for ix, and the vertical axis is for iy, where ix and iy are the indexes of lattice node 
in the x- and y- directions. From top to bottom, the three rows show the density contours at 
times t = 0.5, 0.7, 1, respectively. The common parameters are dx = dy = 0.003, dt = 10~ 5 . 
The collision parameters in Fig. 7(a) are s 5 = s 6 = 10 4 , = s$ = 10 5 , those in Fig. 7(b) 
are s$ = s$ = 10 4 , sj = s$ = 10 4 , and those in Fig. 7(c) are s$ = sq = 10 4 , = s$ = 10 3 , 
10 5 for the others. The collision parameters in Fig. 8(a) are = sq = 10 3 , S7 = sg = 10 5 , 
those in Fig. 8(b) are S5 = = 10 3 , S7 = s$ = 10 4 , and those in Fig. 8(c) are S5 = sq = 10 3 , 
s 7 — s s — 10 3 , 10 5 for the others. Fig. 7 and Fig. 8 show that when the viscosity is 
constant, the small thermal conductivity is beneficial to the development of RM instability. 
Comparing Fig. 7 with Fig. 8, we find when the thermal conductivity is constant, the small 
viscosity is beneficial to the development of RM instability. The thermal conductivity and 
viscosity have inhibition effects on the development of RM instability. Both the unmodified 
model and the modified model get the same results. 



V. CONCLUSIONS 



We propose a MRT Lattice Boltzmann model which works not only for the shocked 
compressible fluids but also for the unshocked compressible fluids. In the new model, a key 
step is the modification of the collision operators of energy flux so that viscous coefficient 
in momentum equation and that in energy equation are consistent no matter if the system 
is shocked or not. The unnecessity of the modification for systems under strong shock is 
analyzed. The new model is validated by some well-known benchmark tests, including (i) 
thermal Couette flow, (ii) Riemann problem, (iii) Richtmyer-Meshkov instability. The first 
system is unshocked and the latter two are shocked. In all the three systems, the Prandtl 
numbers effects are checked. Satisfying agreements are obtained between the new model 



results and analytical ones or other numerical results. Our previous models 
revised in the same way to simulate unshocked compressible flows. 
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